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Abstract 

This paper describes an application of a general purpose computer program, GFSSP 
(Generalized Fluid System Simulation Program) for calculating flow distribution in a 
network of micro-channels. GFSSP employs a f i ni te volume formulation of mass and 
momentum conservation equations in a network consisting of nodes and branches. Mass 
conservation equation is solved for pressures at the nodes while the momentum 
conservation equation is solved at the branches to calculate flowrate. The system of 
equations describing the fluid network is solved by a numerical method that is a 
combination of the Newton-Raphson and successive substitution methods. The numerical 
results have been compared with test data and detailed CFD (computational Fluid 
Dynamics) calculations. The agreement between test data and predictions is satisfactory. 
The discrepancies between the predictions and test data can be attributed to the frictional 
correlation which does not include the effect of surface tension or electro-kinetic effect. 

Nomenclature 


1. Introduction 

Micro-fluidic channel networks [1,2] are being placed on lab-on-a-chip devices, the size 
of a postage stamp, to perform multi-step biochemical processes, or to perform parallel, 
arrayed biochemical experiments. Presently, typical micro fluidic channel networks are 
designed and developed in an iterative process involving very rudimentary modeling, 
followed by several fabrications and testing cycles. This is a costly and ineffective 
approach and often the process does not converge to a suitable system design. The 
reason the iterative process is used is that no commercially available tools exist for 
system-level micro fluidic network design and analysis. The purpose of the present paper 
is to describe the application of a finite-volume based general-purpose computer program 
developed for analyzing propulsion system [3] to model micro-fluidics network. The 
predictions have been compared with experimental data [4] and Navier-Stokes solution 

[5]- 

* 

2. Mathematical Formulation 

The finite volume formulation requires governing equations to be expressed in 
conservative form. The rate of change of a conserved property in a given control volume 
is expressed as the vector sum of transported property from neighboring control volumes 


together with source or sink terms. The unknown variables in the flow circuit of figure 1 
are pressure, temperature and flowrate. These variables are solved from the mass, 
momentum and energy conservation equations in conjunction with correlation of fluid 
friction. 

Figure 1 shows the fluid nodes and branches. There are two types of fluid nodes: 
boundary nodes and internal nodes. The boundary nodes are designated by double 
rectangles (e.g. node 1, 2, 3 etc.). Internal nodes are designated by single rectangles (e.g. 
node 11, 12 etc.). The nodes are connected by branches designated by filled rectangle 
(e.g. 19, 39 etc.). The mass and energy conservation equations are solved at the internal 
nodes and the momentum conservation equations are solved at the branches. 

2.1 Mass Conservation 


The mass conservation equation at the i* node (Figure 2) can be written as 
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Equation 1 implies that the net mass flow from a given node must equate to rate of change 
of mass in the control volume. In the steady state formulation, the left side of the 
equation is zero, such that the total mass flow rate into a node is equal to the total mass 
flow rate out of the node. 


2.2 Momentum Conservation 


The momentum conservation equation at the ij branch can be written as 
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Equation 2 implies that there is a balance between pressure and friction force. The 
frictional force is expressed as a product of friction parameter, Kf and square of flowrate. 


It may be noted that the square of flowrate was expressed as product of m t] and 


m ij 


This 


practice allows the pressure to increase or decrease depending upon the flow direction. It 
may be recalled that the flowrate is a vector property and the pressure is a scalar property. 
The friction parameter, Kf can be expressed as [6] 


K f = 


SJL 
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( 3 ) 


White [9] described a procedure to estimate the friction factor in non-circular duct 
(Figure 3). This procedure consists of the following steps: 
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Figure 1. A typical micro-fluidics network system 


Estimate the hydraulic diameter of the cross-section: 

Dh = (4)(Area) / Perimeter 

Estimate the Poiseuille Number for a particular cross-section. The Poiseuille 
Number can be expressed as a polynomial function of aspect ratio as shown in 
Equation 4. Table 1 provides the coefficients for different geometries. 


Table 1. - Poiseuille Number Coefficients for Non-circular Duct Cross-sections 


Coefficients 

Rectangle 

Ellipse 


Circular 

Section 

Ao 

23.9201 

19.7669 

22.0513 

11.9852 

A, 

-29.436 

-4.53458 

6.44473 

3.01553 

a 2 

30.3872 

-11.5239 

-7.35451 

-1.09712 

a 3 

-10.7128 

22.3709 

2.78999 

0.0 

A4 

0.0 

-10.0874 

0.0 

0.0 


* For b/a < 0.2508 p 0 = ' , where, A 0 = 24.8272, Ai = 0.0479888 

3 . Calculate the friction factor for a non-circular pipe: 

Laminar flow (Re<23001 


4 Po 
Re 


( 5 ) 



Figure 2. Schematic of connections between Nodes by Branches and the indexing 
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Figure 3. Non-circular duct cross-section 


2.3 Energy Conservation Equation 


The energy conservation equation for fluid node i, shown in Figure 2, can be expressed 
following the first law of thermodynamics and using enthalpy as the dependant variable. 
The energy conservation equation based on enthalpy can be written as 
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The MAX operator represents the upwind formulation. 

2.4 Equation of State 

The resident mass in the i th control volume can be expressed from the equation of state 
for a real fluid as 


pV_ 

RTz 


m = 
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For a given pressure and enthalpy the temperature and compressibility factor in equation 
6 is determined from the thermodynamic property program developed by Hendricks et al 
[7]. 

3. Solution Procedure 


The pressure, enthalpy, and resident mass in internal nodes and flowrate in branches are 
calculated by solving equations (1), (6), (7), and (2) respectively. A combination of the 
Newton-Raphson method and the successive substitution method has been used to solve 
the set of equations. The mass conservation (1), momentum conservation (2) and 
resident mass (7) equations are solved by the Newton-Raphson method. The energy 
conservation equations for fluid and solid are solved by the successive substitution 
method. The temperature, density and viscosity are computed from pressure and 
enthalpy using a thermodynamic property program [7], Figure 4 shows the flow diagram 
of the Simultaneous Adjustment with Successive Substitution (SASS) scheme. The 
iterative cycle is terminated when the normalized maximum correction A max is less than 

the convergence criterion C c . A max Is determined from 


A 


max 


= max 
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The convergence criterion is set to 0.001 for all models presented in this paper. The. 
details of the numerical procedure are described in Reference [6]. 


4. Computer Program 

GFSSP (Generalized Fluid System Simulation Program) embodies the mathematical 
formulation and solution procedure described in the previous sections. The program 
structure is shown in Figure 5. The program consists of three modules: Graphical User 
Interface, Solver and User Subroutines. VTASC (Visual Thermofluid dynamics 
Analyzer for Systems & Components) is the Graphical User Interface (GUI). VTASC 
allows user to create a flow circuit using a point and click paradigm. It creates an ASCII 
data file that is read by the solver module and it reads the output data file for post 
processing the results. The solver module reads the data file generated by VTASC. It 
generates all governing equations from network data. The equations are solved by the 
iterative algorithm (SASS). It calls thermodynamic property programs to obtain the 
necessary properties during the iterative cycle. 
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Figure 4. SASS (Simultaneous Adjustment with Successive Substitution) Scheme for 

solving Governing Equations 


5. Computational Fluid Dynamics (CFD) Analysis 

To assist validation of the GFSSP model for micro-fluidics network, a detailed 3-D CFD 
analysis was also carried out. Using CFD, the Navier- Stokes momentum differential 
equations are numerically solved. In this study, the CFD-ACE+ code [5], utilizing finite 
volume numerical method was used for calculating velocity and pressure fields within the 
3-D micro channel configuration. The CFD method utilizes pressure-based SIMPLC 
algorithm and the incompressible flow equations are solved here on a multi-block 
structured grid system. 
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Figure 5. GFSSP Program Structure 
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Figure 6. Micro-fluidic chip considered for model verification 






6. Results & Discussion 


The T-junction micro flow on the Caliper N145 chip was used for validation. Figure 6 
shows the layout and dimensions of micro-fluidic chip considered in the present study. 
The glass chip was micro fabricated by isotropic etching with about 60 micrometer in 
width and 12.4 micrometer in depth. The modeled geometries used for CFD calculations 
are shown in Figure 7 for the cross-section and at the T-junction. The GFSSP model is 
shown in Figure 8. By adjusting the pressure difference among well number 1, 8 and 4, 
different flow fields across the T-junction can be set up. Dyes were introduced in well 
number 1 and well number 8 for flow velocity measurement [4] using vitalization 
technique. To accommodate flow visualization, pressure levels at well number 1, 8 and 4 
were always different with well number 4 being the minimum. For micro-channel with 
this dimension, the electro- viscous effect [1] would be significant for water and low 
concentration ionic solution fluid flows due to the existence of electrical double layer 
(EDL). In addition to pure water, aqueous KC1 solutions of concentration 10' 4 and 10' 2 M 
were used as the testing liquids. Following the studies of [1] and [2], the 10' 2 M 
experiment was taken to be the case without EDL effects when compared to the 
numerical results. 



Figure 7. Modeled cross-section (a) and T-junction (b) geometries for CFD. 

It may be noted that the node numbers in Figure 6 and 8 are not compatible. Node 8 of 
Figure 8 is equivalent to Node 1 of Figure 6. Similarly, Node 11 and 1 of Figure 8 are 
equivalent to Node 8 and 4 of Figure 6. The number of nodes and branches used in the 
numerical model was appropriate for representing resistances in the flow circuit. There 
were two types of resistances in the circuit. The filled rectangle represents flow in a 
rectangular pipe. The micro-channel cross section was assumed rectangular. The other 
type of resistance is denoted by 2K. This option considers the branch as a common 



fittings or valves. The resistance in common fittings and valves can be computed by the 
two-K method [8]. For this option, Ky can be expressed as: 

Tr Ki / Re+ Koo (l + 1 / D) 

K ' (9) 

Where: 

Ki= K for the fitting at Re =1 ; 

Koo = K for the fitting at Re =oo ; 

D = Internal diameter of attached pipe, in. 

The constants Kj and Koo for common fittings and valves are listed in Reference [6]. 



Figure 8. GFSSP model of flow network in chip with micro-channels 
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Figure 9. Comparison between numerical prediction and test data 

Figure 9 shows the comparison between numerical prediction and test data. The velocity 
in branch 34 is plotted with the maximum pressure differential between inlet and outlet. 
Test data has been compared with GFSSP prediction and predictions by CFD-ACE+. It 
may be observed that at low flow velocity, test data, network flow analysis prediction by 
GFSSP and Navier-Stokes flow analysis compare well. However at higher flowrates, 
both predictions deviate from test data. The observed discrepancies could be attributed to 
electro-kinetic effect which has not been accounted for any of the analysis scheme. It 
may also be emphasized that pressure drop is quite accurately predicted by network 
analysis code which is evident from good comparison between GFSSP and CFD-ACE+ 
at all flowrates. 

7. Conclusions 

A numerical method for flow analysis in micro-fluid system has been described. The 
method solves the finite volume formulation of conservation equations for mass, 
momentum and energy in conjunction with friction correlation and thermodynamic 
equation of state. The numerical method is embodied into a general purpose computer 
program for modeling flow distribution in a generic network. The predicted numerical 
solution has been compared with test data as well as multi-dimensional Navier-Stokes 
solution. The network flow analysis results compare well with Navier-Stokes solution at 
different flowrates. Both numerical solutions (network system and Navier-Stokes) 
compare well with test data at low flowrates. Discrepancy with test data exists at higher 
flowrates. The observed discrepancy is attributed to the electro-kinetic effect which has 


not been accounted for. The future research effort is directed to model electro-kinetic 
effect in micro-fluid system. 
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